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ABSTRACT 

This paper focuses on the dynamical implications of close supermassive black hole 
binaries both as an example of resonant phase mixing and as a potential explanation 
of inversions and other anomalous features observed in the luminosity profiles of some 
elliptical galaxies. The presence of a binary comprised of black holes executing nearly 
periodic orbits leads to the possibility of a broad resonant coupling between the black 
holes and various stars in the galaxy. This can result in efficient chaotic phase mixing 
and, in many cases, systematic increases in the energies of stars and their consequent 
transport towards larger radii. Allowing for a supermassive black hole binary with 
plausible parameter values near the center of a spherical, or nearly spherical, galaxy 
characterised initially by a Nuker density profile enables one to reproduce in considerable 
detail the central surface brightness distributions of such galaxies as NGC 3706. 

Subject headings: galaxies: evolution - galaxies: kinematics and dynamics - galaxies: 
structure 
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1. Introduction and Motivation 

Understanding the dynamical implications of a su- 
pcrmassive black hole binary near the center of a 
galaxy is important both because of the insights the 
problem can shed on physical processes associated 
with a time-dependent potential and because, even 
if it itself is not resolvable observationally, the binary 
can have directly observable effects. 

As is well known to nonlinear dynamicists, a time- 
dependent potential can induce significant amounts 
of time-dependent transient chaos, an interval during 
which orbits exhibit an exponentially sensitive depen- 
dence on initial conditions, and resonant couplings be- 
tween the natural frequencies of the time-dependent 
potential and the frequencies of the chaotic orbits 
can trigger efficient resonant phase mixing (Kandrup, 
Vass & Sideris 2003). Like 'ordinary' chaotic phase 
mixing (e.g., Kandrup & Mahon 1994, Merritt & Val- 
luri 1996), this resonant mixing can facilitate a rapid 
shuffling of orbits on different constant energy hyper- 
surfaces. Even more importantly, however, because 
the potential is time-dependent the energies of indi- 
vidual orbits are not conserved, so that resonant mix- 
ing can also facilitate a shuffling of energies between 
different constant energy hypersurfaces. 

For this reason, resonant phase mixing has impor- 
tant implications for collective relaxation in nearly 
collisionless systems (Kandrup 2003), e.g., holding 
forth the prospect of explaining from first principles 
the striking efficacy of violent relaxation (Lynden- 
Bcll 1967) found in simulations and inferred from ob- 
servations (see, e.g., Bertin 2000). That large scale 
collective oscillations could trigger very efficient vio- 
lent relaxation has been shown in the context of one 
simple model, namely orbits of stars in a Plummer 
sphere subjected to a systematic time-dependence 
which eventually damps (Kandrup, Vass & Sideris 
2003). The binary black hole problem provides a 
complementary example of how smaller scale time- 
dependences can also have a surprisingly large effect. 

The binary black hole problem is also interesting 
because the binary can have directly observable con- 
sequences. The fact that energy is not conserved im- 
plies the possibility of readjustments in the density 
profile of stars near the center of a galaxy. In many 
cases this energy nonconservation means that, on the 
average, stars near the center gain energy, which im- 
plies a systematic transport of luminous matter near 
the black holes out to larger radii. To the extent, how- 



ever, that mass traces light, such changes in the den- 
sity distribution translate into predicted changes in 
the observed surface brightness distribution because 
of the presence of such a binary. 

In particular, for reasonable choices of black hole 
masses and orbital parameters, the binary can actu- 
ally cause an 'inversion' in the surface brightness pro- 
file, so that surface brightness is no longer a monoton- 
ically decreasing function of distance from the center. 
Indeed, the simplest models which one might envi- 
sion are adequate to reproduce distinctive features ob- 
served in the brightness distributions of such galaxies 
as NGC 3706, as reported in Laucr et al (2002). 

The first half of this paper, Section 2, considers 
the binary black hole problem as an example of how a 
time-dependent potential can facilitate efficient phase 
mixing in a galaxy. Attention focuses on two sets of 
models, namely the pedagogical example of a constant 
density ellipsoid, corresponding to an anisotropic os- 
cillator potential, and more realistic cuspy density 
profiles consistent with what have been inferred from 
high resolution photometry (e.g., Lauer et al 1995). 

One important issue here involves determining as 
a function of amplitude (i.e., black hole masses) 
and frequency (i.e., orbital period) when the time- 
dependent perturbation can have a significant effect. 
A second involves determining the degree to which the 
efficacy of energy and mass transport reflect the de- 
gree of chaos exhibited by the orbits, both in the pres- 
ence and the absence of the perturbation. To what ex- 
tent, e.g., does efficient energy transport require that 
a large fraction of the orbits in the time-dependent 
potential be chaotic? Does resonant phase mixing 
rely crucially on the presence of transient chaos? 

Another issue involves determining the extent to 
which the bulk manifestations of a black hole binary 
vary for spherical, axisymmetric, and nonaxisymmet- 
ric (e.g., triaxial) galaxies. Is it, e.g., true that spher- 
ical and nearly spherical systems are impacted less 
by the presence of a supermassive binary since, in the 
absence of the binary, all or almost all of the orbits 
are regular? In a similar vein, one would like to un- 
derstand the extent to which the effects of the binary 
depend on the steepness of the cusp. And, perhaps 
most importantly, it would seem crucial to determine 
how the size of the 'sphere of influence' of the binary 
depends on the size of the black hole orbits and their 
masses. Perhaps the most important conclusion here 
is that this 'sphere' can be much larger than the size 
of the black hole orbits. For plausible choices of pa- 
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rameter values, black holes moving along orbits with 
size ~ Th can significantly impact the density distri- 
bution at radii as large as ~ 10 — 20r/j or more. 

All these issues have important implications for 
determining when a supermassive black hole binary 
might be expected to have observable consequences. 
The second half of the paper, Sections 3 and 4, con- 
siders these consequences. Section 3 considers the 
generality of the simple models considered in Section 
2, which assume circular orbits and equal mass black 
holes, and then focuses on direction-dependent effects 
which must be understood to determine how poten- 
tially observable quantities depend on the relative ori- 
entation of the observer and the binary. 

Section 4 focuses in detail on one specific observ- 
able prediction, namely that supermassive black hole 
binaries can alter the density distribution near the 
center of a galaxy. This involved: (i) generating TV- 
body realisations of density distributions consistent 
with a Nukcr Law (Laucr et al 1995); (ii) evolving 
these iV-body systems in the fixed time-dependent 
potential corresponding to the galaxy plus orbiting 
black holes; (iii) determining how the initial density 
distribution changes over the course of time; and, (iv) 
presuming that mass traces light, integrating the re- 
sulting density distribution along the line of sight to 
obtain a surface brightness profile. These are not self- 
consistent computations; but they can at least provide 
strong indications as to what the expected effects of 
the binary would be. The crucial point, then, is that 
such an exercise results generically in brightness dis- 
tributions that resemble qualitatively the forms re- 
ported in Lauer et al (2002); and that by fine-tuning 
parameters within a reasonable range, one can repro- 
duce many of the details of what is actually observed. 

Section 5 summarises the principal conclusions and 
discusses potential implications. 

2. Dynamical Effects of Supermassive Black 
Hole Binaries 

2.1. Description of the experiments 

The computations described here involved orbits 
evolved in potentials of the form 



V(x,y,z) = V (x,y,z)- 



M 



M 



|r-n(i)| |r-r 2 (t)| 



. (1) 



where Vo is time-independent and ri and r 2 corre- 
spond to circular orbits in the x — y plane, i.e., 

x\{t) — rh sinwi, Vi(t) = Th coswt, z\(t) = 0, 

(2) 

and r 2 = — ri. Some of the computations focused on 
a harmonic oscillator potential, 



Vo(x,y, z) = -m 2 , 



(3) 



with 



m 2 = [x/af + (y/bf + {z/cf. (4) 
Others focused on more realistic potentials of the form 
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with 7 is a cusp index, assumed to satisfy < 7 < 2. 
The axis ratios a, b, and c were selected of order unity. 

The assumptions that the black holes are in circu- 
lar orbits and that they have equal masses might ap- 
pear an extreme idealisation. However, as will be dis- 
cussed in Section 3, it appears that relaxing these as- 
sumptions does not change the principal conclusions. 
This model appears structurally stable towards modest 
changes in the orbital parameters of the binary. 

For a = b = c = 1, cq. (5) reduces to the spherical 
Dchnen (1994) potential with unit mass, and, quite 
generally, for large r, V — ► — 1/m. Thus for axis ra- 
tios of order unity, one can interpret eq. (5) as the 
potential for a galaxy with mass M g rts 1.0. For non- 
spherical systems, eq. (5) yields density distributions 
different from Merritt and Fridman's (1996) triaxial 
Dehnen models, in that it is V, rather than p, that is 
constrained to manifest ellipsoidal symmetry. 

This potential is unrealistic in that, for large radii, 
V does not become spherically symmetric; and one 
can also argue that it is unrealistic in the sense that, 
assuming mass traces light, the isophotes become 
peanuty for axis ratios far from spherical. Given, 
however, that one is interested primarily in physical 
processes in the central portions of the galaxy, the 
r — > 00 asymptotic behaviour is largely unimportant; 
and it should be recalled that the isophotes in 'real' 
galaxies tend to manifest systematic deviations from 
ellipticity (e.g., Kormendy & Bender 1996). This po- 
tential has the huge advantage that, unlike Merritt 
and Fridman's nonspherical Dchnen potential, it can 
be expressed analytically, thus reducing by two orders 
of magnitude or more the time required for orbital 
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integrations. Moreover, as discussed in Section 4, for 
the case of spherical symmetry the behaviour of or- 
bits in this potential is very similar to orbits evolved 
in the potential associated with a Nuker Law, at least 
for those choices of Nuker parameters for which the 
potential can be expressed analytically. 
For spherical systems with a = b = c = 1 , 

M(r) = [r/(l + r)] (3 ~ 7) (6) 

is the mass within r of the galactic center. For axis 
ratios of order unity, eq. (7) also provides a reasonable 
estimate for moderately nonspherical systems. 

The computations here involved black hole masses 
in the range 0.005 < M < 0.05 and radii satisfying 
0.005 < r h < 0.5. Following Merritt and Fridman's 
(1996) normalisation for their triaxial Dehncn model, 
one can translate the dimcnsionless model into phys- 
ical units by defining the correspondence 

t = l 1.46xl0 6 M 1 - 1 1/2 ajl£yr, (7) 

where 

Mn = (io^k) and flkpc = (ife) • 

One can identify an energy-dependent dynamical 
time to following either Merritt and Fridman, who 
related it to the period of a specific type of regular or- 
bit, or Kandrup & Siopis (2003), who proposed a pre- 
scription based on the times between turning points 
in representative orbits. Those two prescriptions yield 
results in agreement at the 5% level or better. More 
generally, for axis ratios of order unity, at least for 
small radii the angle-averaged density and energy dis- 
tributions are relatively similar to those associated 
with 'true' maximally triaxial Dchnen models, so that 
a dynamical time to(E) can be estimated to within 
20% or so from Table 1 in Merritt & Fridman (1996). 
This implies, e.g., that, for 7 = 1.0, a time t = 512 
corresponds to roughly lOOt^ for stars in the 20% 
mass shell or, equivalently, ~8x 10 8 M 11 1 ^ 2 a^p^ yr. 

A 'realistic' value for the frequency can be esti- 
mated easily. Suppose, e.g., that a = b = c = 1.0. If 
M <C M(rh), with M(r/ l ) the galactic mass contained 
within radius rh, the black holes can be viewed as test 
particles moving in the galactic potential. This im- 
plies that 

^ 2 -^ 7 (l + ^) 7 " 3 ; (8) 

so that, e.g., uj = r h 1 ^ 2 (l+rf l )~ 1 for 7 = 1.0. If, alter- 
natively, M » M(rh) the potential associated with 



the galaxy can be neglected and one is reduced de 
facto to the circular, equal mass two-body problem, 
for which uj = ^jMj\r\. For 7 = 1.0, M = 0.01, and 
rh = 0.05, fiducial values considered in many of the 
computations, the galactic potential can be neglected 
in a first approximation, so that w w >/20 w 4.47. 

However, for much of this section, w was viewed 
as a free parameter so that, for fixed amplitude and 
geometry, one can explore the response as a function 
of driving frequency. This enables one to determine 
the extent to which the response manifests a sensitive 
dependence on frequency, which can provide impor- 
tant insights into the resonant couplings generating 
the response. 

Attention focused primarily on the statistical prop- 
erties of representative orbit ensembles, integrated 
from sets of > 1600 initial conditions. These were 
generated by uniformly sampling a specified constant 
energy hypersurface as defined in the limit M — > 
using an algorithm described in Kandrup & Siopis 
(2003). Allowing for the black holes changes the 
initial energies, so that one is de facto sampling 
a 'slightly thickened' constant energy hypersurface. 
The initial conditions were integrated forward for a 
time t = 512. The integrations also tracked the 
evolution of a small initial perturbation, periodically 
renormalised in the usual fashion {e.g., Lichtenberg 
& Lieberman 1992), so as to extract estimates of the 
largest finite time Lyapunov exponent. 

For each simulation, specified by a, b, c, M, rh, 
and u, following quantities were extracted: 

(i) the fraction / of 'strongly chaotic' orbits, esti- 
mated as in Kandrup & Siopis (2003) (as discussed 
in Kandrup, Vass & Sideris [2003], because the po- 
tential is time-dependent it is often difficult to make 
an absolute distinction between regular and chaotic 
orbits, although it is relatively easy to identify orbits 
that are 'strongly chaotic'); 

(ii) the mean value (\) of the finite time Lyapunov 
exponents for the strongly chaotic orbits; 

(iii) the mean value (SE) of the energy shift SE = E{t) 
—E(0) for all the orbits at various times t > 0; and 

(iv) the dispersion o-se associated with these shifts. 
The data were also analysed to determine the func- 
tional forms of (SE(t)) and crsE(t), and to search for 
correlations between changes in energy and values of 
finite time Lyapunov exponents for individual orbits 
within a single ensemble. 

Other integrations tracked phase mixing in initially 
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localised ensembles, so as to determine the extent to 
which such mixing resembles ordinary chaotic phase 
mixing in a time-independent potential (e.g., Mer- 
ritt & Valluri 1996, Kandrup 1998) or resonant phase 
mixing in a galaxy subjected to large scale bulk oscil- 
lations (Kandrup, Vass & Sideris 2003). 

2.2. Statistical properties of orbit ensembles 

Overall, as probed by the shuffling of orbital en- 
ergies, there is a broad and comparatively efficient 
resonant response. For fixed values of a, b, c, M, and 
r/j, the range of 'interesting' frequencies lo can be two 
orders of magnitude or more in breadth. One does not 
need to 'fine-tune' lo to trigger an efficient shuffling of 
energies. However, the resonance can exhibit substan- 
tial structure, especially for the case of a spherically 
symmetric oscillator potential. Superimposed upon a 
smooth overall trend, quantities like (SE) can exhibit 
a complex, rapidly varying dependence on lo. 

Consider, e.g., Fig. 1, which plots (SE) and o~se as 
functions of lo at time t — 512 for two oscillator mod- 
els, one spherical and the other triaxial. Both have 
M = 0.05 and r/j = 0.3 The curves for the two mod- 
els have envelopes with a comparatively simple shape 
but, for the spherical model, an enormous amount 
of substructure is superimposed. This substructure 
reflects the fact that all unperturbed orbits oscillate 
with the same frequency. Indeed, close examination 
reveals that the resonances are associated with in- 
teger and (to a lesser degree) half-integer values of 
lo, harmonics of the unperturbed natural frequency 
lo = 1.0. In an axisymmetric system, there are two 
natural frequencies, which can yield a yet more com- 
plex response pattern. If, however, M and r^ are 
chosen large enough to elicit a significant response, 
the resonances typically broaden to the extent that 
much, if not all, that structure is lost. Allowing for a 
fully triaxial system leads to three unequal frequen- 
cies, which yields such a plethora of harmonics that, 
even for comparatively weak responses, the short scale 
structure is largely lost. 

It is evident from Fig. 1 that, although (SE) is sub- 
stantially larger for the triaxial than for the spherical 
model, asE is comparable. What this means is that, 
even though the spherical model leads to a smaller 
systematic shifting in energies, the energies of orbits 
in these two models are shuffled to a comparable de- 
gree. The observed differences in (SE) do not re- 
flect the fact that the nonsphcrical model is triaxial. 
Rather, they appear again to reflect the fact that, for 



a spherical system, there is only one characteristic fre- 
quency for the unperturbed orbits. Modest deviations 
from spherical symmetry, be these either axisymmet- 
ric or not, suffice typically to yield amplitudes more 
closely resembling Fig. 1 (c) than 1 (d). 

As indicated in Fig. 2, there is a threshold value of 
M below which no substantial response is observed; 
and, similarly, the response 'turns off for higher- 
energy orbits that spend most of their time far from 
the binary. However, for 'interesting' choices of M, 
as probed, e.g., by (SE(lo)) or (tse(w), the resonance 
has a characteristic shape. As the frequency increases 
from lo = 0, (SE) and ose exhibit a rapid initial in- 
crease, peak at a maximum value, and then begin a 
much slower decrease. For fixed parameter values, the 
value of the frequency triggering the largest response 
is roughly independent of mass, but it is true that, 
for larger M, the relative decrease in (5E(lo)) and 
o~se(u) with increasing lo is slower than for smaller 
M. This is consistent with the notion that, for larger- 
amplitude perturbations, higher-order harmonics be- 
come progressively more important. Fig. 2 also shows 
the peak frequency is a decreasing function of rv In 
particular, when the black holes are closer together 
one requires larger lo to elicit a significant response. 

Efficient shuffling of energies seems tied unambigu- 
ously to the presence of large amounts of chaos, as 
probed by the fraction / of strongly chaotic orbits 
and, especially, the size of a typical finite Lyapunov 
exponent (\) ■ Large / and (x) do not guarantee large 
changes in energies, but they are an essential prerequi- 
site. In some cases, notably nearly spherical systems, 
/ and (x) are very small in the limit lo — > 0. As w 
increases, however, / and especially (x) also increase; 
and, for values of lo sufficiently large to trigger an 
efficient response, the ensemble will be very chaotic 
overall. For values of u> in the resonant region, / and 
(x) tend to exhibit only a comparatively weak depen- 
dence on lo. 

Orbit ensembles evolved in spherical, axisymmet- 
ric, and triaxial Dehnenesque potentials exhibit res- 
onance patterns quite similar both to one another 
and to the patterns observed in nonspherical oscil- 
lator models, although some relatively minor differ- 
ences do exist. Fig. 3 exhibits data for two such mod- 
els, one spherical and the other triaxial, each with 
7 = 1.0, rh = 0.05, and M = 0.01. The models were 
both generated for ensembles of initial conditions with 
E = -0.70 and (r in ) w 0.33. The black hole radius 
r/j = 0.05 corresponds roughly to the 0.2% mass shell. 
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The triaxial oscillator model in Fig. 1 and the 
Dehnenesque models in Fig. 3 are representative in 
that modest changes in the parameters of the binary 
do not lead to significant qualitative changes in the 
response. Equally important, however, they are also 
robust towards changes in axis ratio. Axisymmet- 
ric and slightly triaxial Dehnenesque models (e.g., as 
nonspherical as a 2 = 1.05, b 2 = 1.00, and c 2 = 0.95) 
yield results very similar to the spherical case. Other 
'strongly' triaxial models yield results similar to the 
particular triaxial model exhibited here. 

Even for spherical and nearly spherical systems, 
the relative measure / of strongly chaotic orbits tends 
to be large even for u = 0, this corresponding to 
stationary but separated black holes. One does not 
require a strong time- dependence to generate a large 
measure of chaotic orbits. However, it is true that, for 
axisymmetric and other nearly spherical systems, the 
size of a typical Lyapunov exponent tends to be con- 
siderably smaller than for strongly triaxial systems. 

In significantly triaxial models, the degree of chaos, 
as probed by / or (x) , is a comparatively flat function 
of ui. By contrast, in nearly spherical and axisymmet- 
ric systems, the degree of chaos, especially as probed 
by (x), increases rapidly with increasing u) until it 
becomes comparable to the degree of chaos exhibited 
by strongly triaxial models. The obvious inference 
is that, when the system is nearly spherical or ax- 
isymmetric, the time- dependence associated with the 
orbiting binary is required to give the chaotic orbits 
particularly large Lyapunov exponents. 

An analogous result holds for the mean energy shift 
(SE). Quite generally, (SE) — > for u) — > and in- 
creases with increasing u. However, the initial rate of 
increase is typically much larger for significantly triax- 
ial models than for nearly spherical and axisymmetric 
systems. This has an important practical implication: 
Because larger frequencies are required to trigger the 
resonance in galaxies that are nearly axisymmetric, in 
nearly spherical or axisymmetric galaxies black holes 
of given mass must be in a tighter orbit before they 
can trigger a significant response. 

2.3. Shuffling of energies as a diffusion pro- 
cess 

Overall, the shuffling of energies induced by the 
black hole binary is diffusive, although the basic pic- 
ture depends on the amplitude of the response. When 
changes in energy experienced by individual orbits are 



relatively small, the dispersion tends to grow diffu- 
sively, i.e., asE oc t 1 / 2 . The mean shift in energy typ- 
ically grows more quickly, being reasonably well fit by 
a linear growth law (SE) oc t. Alternatively, when the 
response is stronger, it is the mean shift that grows 
diffusively, i.e., (SE) oc t 1 / 2 , whereas the dispersion is 
well fit by a growth law ctse c< t 1 / 4 . Examples of both 
sorts of behaviour can be seen in Fig. 4. 

One might have supposed that, since the shuffling 
of energies is associated with the presence of chaos, 
changes in energy should grow exponentially. This 
however, does not appear to be the case, at least 
macroscopically. The initial response of the orbits 
(t < 5 to or so) may be exponential, but it is evi- 
dent that, overall, the response is diffusive. Time- 
dependent chaos does not trigger exponentially fast 
mixing in energies. However, it can still be extremely 
important in that it allows comparatively efficient 
shufflings of energies that would be completely impos- 
sible in a time-independent Hamiltonian system. 

One final point should be stressed. That changes in 
energy are diffusive, reflecting a slow accumulation of 
energy shifts, corroborates a fact also evident from an 
examination of individual orbits: Changes in energy 
experienced by individual orbits do not result from 
single close encounters with the black holes. Instead, 
they really do reflect resonance effects associated with 
the time-dependent potential. 

2.4. Correlations amongst orbital properties 
for different orbits within an ensemble 

Orbits with smaller finite time Lyapunov expo- 
nents x tend to exhibit energy shifts that are smaller 
in magnitude \SE\. Orbits with large x can experi- 
ence both large and small net changes in energy. As 
has been observed in other time-dependent potentials 
(Kandrup & Terzic 2003), the fact that an orbit is 
chaotic does not necessarily imply that it will exhibit 
large, systematic drifts in energy over a finite time in- 
terval. However, the energy shifts in orbits with small 
X are invariably small. 

When the response is weak, so that the dispersion 
of the ensemble evolves diffusively, changes in energy 
exhibited by individual orbits are comparably likely to 
be positive or negative. However, when the response 
is stronger, so that the mean shift evolves diffusively, 
the energies of individual orbits tend instead to in- 
crease systematically. 

When the response is relatively weak and changes 
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in energy are equally likely to be cither positive or 
negative, the distribution of energy shifts n(SE) is 
typically well fit by a Gaussian with mean roughly 
equal to zero. However, when the response becomes 
stronger, n(SE) becomes distinctly asymmetric and 
cannot be well fit by a Gaussian, even allowing for a 
nonzero mean. 

Correlations between the 'degree' of chaos and the 
'degree' of energy shuffling experienced by individ- 
ual orbits are perhaps best illustrated by extract- 
ing energy shifts 6E at different times ti for indi- 
vidual orbits, computing the mean and dispersion, 
(SE) and <jse, associated with the resulting time se- 
ries {SE(ti)}, and demonstrating how these moments 
correlate with the value of the finite time Lyapunov 
exponent \- Examples of such an analysis are exhib- 
ited in Fig. 5. The obvious point is that the moments 
are invariably small when \ is small, whereas larger 
X typically implies larger values of \{SE}\ and use- 

2.5. Chaotic and resonant phase mixing 

The time-dependent potential associated with the 
black hole binary can alter 'ordinary' chaotic phase 
mixing in at least two important ways. 

The time-dependent potential tends to increase 
both the fraction of chaotic orbits and the size of 
a typical Lyapunov exponent. If a galaxy is in a 
(nearly) time-independent (near- Equilibrium state, 
the relative measure of (at least strongly) chaotic or- 
bits should be relatively small, since presumably one 
requires large measures of regular (or nearly regu- 
lar) orbits to provide the 'skeleton' of the interest- 
ing structures associated with those chaotic orbits 
which are present (Binney 1978). Introducing a time- 
dependent perturbation leads oftentimes to a signifi- 
cant increase in the relative measure of chaotic orbits. 
Moreover, even when the time-dependence does not 
significantly increase the measure of chaotic orbits, it 
can make already chaotic orbits more unstable, thus 
allowing them to mix more efficiently. 

Because energy is no longer conserved, the time- 
dependent potential also allows mixing between dif- 
ferent constant energy hypersurfaces, which is com- 
pletely impossible in the absence of a time-dependence. 
Overall, this mixing of energies is not as efficient a 
process as mixing in configuration or velocity space. 
However, the resonant mixing of energies associated 
with chaotic orbits still plays an important role. 

An example of such resonant phase mixing is pro- 



vided in the left panels of Fig. 6, which track an ini- 
tially localised ensemble with E = —0.70 in a spher- 
ical Dehnen potential with 7 = 1 and a = b = c = 

1.0, allowing for black hole parameters M — 0.005, 
rh = 0.05, and ui — -\/l0. The right panels track the 
same ensemble, evolved identically except that uj = 0. 
Two points are immediate. One is that, for the re- 
alistic case when u> ^ 0, a time t — 64, correspond- 
ing to <~ 10 8 M 1 ~ 1 1 ^ 2 a||p^ yr, is sufficient to achieve a 
comparatively well mixed configuration. Achieving a 
comparable degree of mixing for the u = system 
requires a time t > 512. The other point is that or- 
bits in the ensemble evolved with have diffused 
to radii r > 0.3, which is impossible for orbits in the 
uo = ensemble, for which energy is conserved. 

3. Observational Consequences of the Dy- 
namics 

3.1. Generality of the idealised model 

Attention hitherto has focused on the dynamical 
consequences of a supermassive black hole binary, 
viewed as the prototype of a time-dependent pertur- 
bation acting in a galaxy idealised otherwise as a col- 
lisionless equilibrium. The object of this and the fol- 
lowing section is to consider instead potentially ob- 
servable consequences, the most obvious of which is a 
changing surface brightness distribution induced by a 
readjustment in the mass density as stars are trans- 
ported to larger radii. 

In so doing, one can proceed by viewing the host 
galaxy as a superposition of orbit ensembles with dif- 
ferent energies E and, for various choices of binary 
parameters, determining when, for any given value of 
E, the binary can have an appreciable effect, e.g., by 
generating a large energy shift (SE) . As described al- 
ready, the response will only be large when the size r^ 
of the binary orbit is sufficiently small that the total 
black hole mass Mi + M 2 > M{rh)- This, however, 
implies that, in a first approximation, the frequency 
of the binary can be estimated neglecting the bulk po- 
tential of the galaxy. Thus, relaxing the assumptions 
of equal masses and strictly circular orbits, 

/ M 1 +M 2 
U = V ^ > ( 9 ) 

with A the value of the semi-major axis. 

Perhaps the most obvious question here is simply: 
For fixed E and A, how do quantities like (SE) de- 
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pcnd on the total mass M to t — Mi + M 2 ? The an- 
swer is that, at least for 'realistic' binary black hole 
masses, i.e., Mi and M 2 < 0.01M ga ;, (SE) is a com- 
paratively smooth, monotonically increasing function 
of M to t- For very small masses, there is essentially 
no response; but, beyond a critical mass, the pre- 
cise value of which depends on properties of the host 
galaxy, the dependence is roughly power law in form, 
i.e., (SE) oc Mf ot , with the power p typically in the 
range 1 < p < 2. Examples of this behaviour are 
exhibited in the left panels of Fig. 7, which show the 
effects of increasing the total mass for five different 
models, one spherical, one prolate axisymmetric, one 
oblate axisymmetric, and two genuinely triaxial. This 
particular set of examples again incorporated circular 
orbits and equal black hole masses; but, as will be 
discussed below, these assumptions are not crucial. 

A second obvious question is: How small must 
the binary orbit be in order to elicit a significant re- 
sponse? Physically, one might suppose that the bi- 
nary was initialised in a comparatively large orbit as 
the result of a merger of two colliding galaxies; but 
that the orbit slowly decayed via dynamical friction, 
allowing the black holes to sink toward the center of 
the galaxy. However, within the context of such a 
scenario the crucial issues to determine are (i) when 
the binary can begin to have a large effect, i.e., how 
small the orbit must be; and (ii) when the effects of 
the binary 'turn off' again. These issues are addressed 
in the right panels of Fig. 7, which exhibit (SE) as a 
function of ru for the same five galactic models used 
to generate the left panels. 

Two points are evident: (1) The binary has its 
largest effect when ru is substantially smaller than the 
typical radius of the orbits with the specified energy. 
The ensembles considered were each comprised of or- 
bits with initial energy E = —0.70 and mean radius 
(r) i=a 0.33, but the maximum response was observed 
for r/j ~ 0.04, i.e., a size roughly ten times smaller! 
This reflects the fact that mass and energy transport 
have been triggered by a resonance, rather than by 
direct binary scatterings of individual stars with the 
black holes. One needs a very tight binary orbit to 
get frequencies sufficiently large to trigger a signifi- 
cant response. (2) As noted already, for the triaxial 
models the effects of the binary 'turn on' at substan- 
tally larger values of Th than for the spherical and ax- 
isymmetric systems. This would suggest that a black 
hole binary could have an especially large effect in a 
strongly triaxial galaxy: since the range of black hole 



sizes that can have an appreciable effect is substan- 
tially larger, the time during which the resonance will 
act should presumably be longer. 

But how generic are the idealised computations de- 
scribed in Section 3? It might not seem unreasonable 
to assume that the black holes follow nearly circular 
orbits, since dynamical friction will tend to circularise 
initially eccentric orbits; but the assumption of equal 
mass black holes is clearly suspect. 

Computations show that varying the eccentricity 
e within reasonable bounds has only a comparatively 
minimal effect. Increasing e from values near zero to 
a value as large as e = 0.5 will not change quantities 
like (SE) by more than 25%; and, in general the effect 
is much smaller even than this. This is, e.g., evident 
from the left panels of Fig. 8, which were generated 
for the same five models considered in Fig. 7. 

As is evident from the right hand panels of Fig. 8, 
there is a substantially stronger, systematic depen- 
dence on the mass ratio M2/M1. For fixed M tot = 
Mi + M 2 , the largest effects arise for Mi w M 2 ; but 
even here the dependence on the mass ratio is not 
all that sensitive. In particular, for all but the tri- 
axial models, the response is a relatively flat func- 
tion of M 2 /M to t for M 2 /M tot ^0.25, so that, for fixed 
Mi + M 2 , mass ratios 1/3 £ M 2 /M x £ 1 yield com- 
parable results. It is true that, for fixed semi-major 
axis A and total mass M tot , the effect of the binary 
is significantly reduced for M 2 <C Mi, but the rea- 
son for this is obvious: When M 2 <C Mi, the more 
massive black hole is located very near the center of 
the galaxy. This implies, however, that, even if the 
binary has a very high frequency, the more massive 
black hole remains too close to the center to have an 
appreciable effect at large radii. The smaller black 
hole is typically found at much larger values of r, but 
its mass is too small to have a significant effect. 

3.2. Systematic changes in density 

Changes in energy induced by transient chaos lead 
generically to a readjustment in bulk properties like 
density; and, to the extent that there is an average 
increase in energy, this readjustment implies a sys- 
tematic displacement of stars to larger radii. To see 
how this effect can proceed, one can sample a constant 
energy hypersurface to generate a set of initial con- 
ditions, integrate those initial conditions into the fu- 
ture, and then compare angle-averaged radial density 
distributions p(r) generated at various times t > 0. 
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The left panels of Fig. 9 summarise results for a 
model with a 2 = 1.25, b 2 = 1.0, and c 2 = 0.75, 
assuming circular orbits with Mi — M 2 = 0.01, 
Th = 0.05, and oj = \/20. The ensemble was so con- 
structed that E = -0.70 and (r m ) w 0.33. The five 
panels exhibit the density distributions at t — 0, 16, 
32, 64, and 128, the last corresponding physically to 
~ 2x 10 8 M n 1 ^ 2 a^ yr. The right panels exhibit anal- 
ogous data for the same ensemble and potential but 
with the black holes held fixed in space, i.e., oj = 0. 

The density distribution remains essentially un- 
changed for the time-independent uj = potential, 
but the realistic case with u = >/20 leads to a signifi- 
cant density readjustment. (Minor changes in the lu = 
model reflect a modest readjustment to the insertion 
of the fixed binary in an equilibrium generated with- 
out a binary.) Already by t = 16 (11 binary periods), 
corresponding to an interval <~ 2.5x 10 7 M 1 ^ L 1 ^ 2 a^ yr, 
there is a pronounced decrease in density in the range 
0.3 < r < 0.5 and an increase in density at larger 
radii. Initially the trajectories are restricted energet- 
ically to r < 0.6. By t = 128, more than 13% of the 
trajectories are located at r > 1.0. 

3.3. The size of the 'sphere of influence' 

Figure 9 demonstrates that a black hole binary can 
significantly impact orbits which spend most of their 
times at radii 3> rv The obvious question, however, 
is: how much larger? To answer this question one 
can evolve ensembles with a variety of different ini- 
tial radii, and determine their response as a func- 
tion of r. The results of two such investigations are 
summarised in Fig. 10. In each case, the configura- 
tion corresponded to a spherical Dehncn model with 
a = b = c = 1.0 and a binary with M 1 = M 2 = 0.005, 
r/j = 0.25, and u> = 0.2828. The left panels are for a 
model with 7 = 0.0; the right panels for 7 = 1.0. 

The 'sphere of influence' is in fact quite large, ex- 
tending out to r > 4, even though ru = 0.25. More- 
over, it is evident that the ensembles which expe- 
rience the most shuffling in energies, as probed by 
and o~se, are precisely those ensembles with 
the largest Lyapunov exponent (x). Indeed, for the 
7 = 0.0 and 7 = 1.0 models, the rank correlation 
between the mean shift and the mean expo- 

nent (x) for different ensembles, are, respectively, 
H{{5E), (x» = 0.615 and 0.613. 

It is also clear that the value of the cusp index 7 
has a significant effect on the details of the response. 



The value of 7 does not have a large effect on the size 
of the binary 'sphere of influence', but it does impact 
the amplitude of the response and how that response 
correlates with radius. In both cases, there is a signif- 
icant response for 0.15 < r < 6.0, but the response in 
this range, as probed both by the degree of shuffling 
in energies, is somewhat larger for the cuspy model. 
Even more strikingly, however, the cusp appears to 
reduce both the size of the Lyapunov exponents and 
the degree of shuffling at very small radii. For the 
cuspy model with 7 = 1.0, comparatively little shuf- 
fling of energies and comparatively small amounts of 
chaos are observed at radii <C r^. The very lowest 
energy stars tend to be more regular and to be less 
susceptible to resonant mixing. 

3.4. Anisotropy 

To what extent does the mass transport induced 
by a supermassive black hole binary depend on di- 
rection? Even if, e.g., the host galaxy is modeled 
as exactly spherical, the binary breaks the symmetry 
and, as such, could introduce anisotropies into a com- 
pletely isotropic ensemble of stars. This is important 
since such anisotropies would imply that changes in 
visual appearance induced by the binary could depend 
appreciably on the observer's viewing angle. 

As a simple example, one can consider the direction- 
dependent density distributions associated with a 
uniform sampling of a constant energy hypcrsur- 
face which, assuming a spherical potential, implies 
a spherically symmetric density distribution and an 
isotropic distribution of velocities. One example 
thereof is exhibited in Fig. 11, which was generated 
for a 7 = 1.0 Dehnen model with a — b = c = 1 con- 
taining a binary executing a circular orbit in the x — y 
plane with the 'correct' Kepler frequency. Here the 
top left panel exhibits spatial distributions at times 
t = and t = 512; the top right panel shows the 
corresponding velocity distributions. At t = 0, the 
spatial and velocity distributions arc equal modulo 
statistical uncertainties; at late times they differ sys- 
tematically, but it remains true that n(|a;|) = «(|y|) 
and n(|ua;|) = nduyl). There is clearly a systematic 
outward transport of stars in all three directions, but 
it is also evident that there is a larger net effect on 
the spatial components in the plane of the orbit. Simi- 
larly, there is a modest shift in velocities which, again, 
is more pronounced in the x and y components. The 
bottom two panels contain plots of, respectively, the x 
and y and the x and z coordinates at t = 512. These 
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panels confirm that the final distribution is more ex- 
tended in the plane of the binary than in the orthog- 
onal direction. In particular, it could easily be misin- 
terpreted as a disc or a torus. 

But what if the host galaxy is already nonspher- 
ical? If, e.g., the galaxy is genuinely triaxial, one 
might suppose that the binary will have settled into a 
symmetry plane; but, assuming that this be the case, 
there are at least two obvious questions. (1) How 
does the overall response depend on which symmetry 
plane? (2) For a binary oriented in a given plane, 
to what extent do observable properties depend on 
viewing angles? 

Both these questions were addressed as before by 
evolving uniform samplings of constant energy hyper- 
surfaces which yield a triaxial number density but are 
still characterised by an isotropic distribution of ve- 
locities. The result of one such computation is sum- 
marised in Fig. 12. Here panel (a) exhibits the initial 
density distributions; the remaining three panels ex- 
hibit the corresponding distributions at t = 512 for 
three different integrations, with the binary oriented 
in the x — y, y — z, and z — x planes. 

Overall the 'angle-averaged' properties of the dif- 
ferent simulations are very similar: The mean short 
time Lyapunov exponents (x) for the three different 
runs agree to within 10%, and even smaller variations 
were observed for quantities like (SE). Indeed, the 
shape of the galaxy seems more important than the 
orientation of the binary. For all three binary orien- 
tations, one observes that the largest effect is in the 
x-direction, which corresponds to the long axis, and 
the smallest in the short-axis z-direction. The details 
of the response observed here depend to a consider- 
able extent on both the shape of the potential and the 
energy of the initial ensemble. In particular, for some 
choices the response is largest in the short-axis rather 
than long-axis direction. However, it seems true quite 
generally that the orientation of the binary is compar- 
atively unimportant. There remains a dependence on 
viewing angle but, if anything, this effect is somewhat 
weaker than for the case of spherical systems. 

For axisymmetric systems with the binary oriented 
in the x—y symmetry plane, one finds generically that 
distributions in the x and y directions agree to within 
statistical uncertainties, but that the z-direction dis- 
tributions differ systematically. In some cases (de- 
pending on both shape and binary parameters), there 
is more mass transport in the z direction; in others the 
effect is more pronounced in the x and y directions. 



These differences likely reflect the fact that this mass 
transport is triggered by a resonance. The unper- 
turbed orbits have different characteristic frequencies 
in different directions, but this would suggest that the 
resonant coupling could well be stronger (or weaker) 
in one direction than in another. 

4. Modeling Luminosity Profiles in Real Gal- 
axies 

4.1. Basic strategy 

The objective here is to show that the physical ef- 
fects discussed above, seemingly the inevitable con- 
sequence of a supermassive black hole binary in the 
center of a galaxy, could provide a natural explana- 
tion of the fact that, in a number of galaxies that 
have been observed using WFPC 2 {e.g., Lauer et al 
2002), the projected surface brightness distribution 
in a given direction is not a monotonically decreasing 
function of distance from the center of the galaxy. 

The computations described here are not com- 
pletely realistic. As in Section 2, they assume black 
holes of exactly equal mass executing exactly circu- 
lar orbits, and the computed orbits of test stars are 
not fully self-consistent since one is neglecting both 
changes in the form of the bulk potential as stars are 
displaced from their original trajectories and the slow 
decay of the binary orbit. They do, however, demon- 
strate that allowing for a binary of relatively small 
size, ~ 10 pc, comprised of black holes with mass 
<; 1% the mass of the galaxy, leads generically to lu- 
minosity dips of the form that have been observed. 
Moreover, fine-tuning parameters within a reasonable 
range of values allows for the possibility of a compar- 
atively detailed (albeit in general nonunique) fit to 
observations of specific galaxies. 
The basic programme is as follows: 

• Generate iV-body realisations of a spherical galaxy 
characterised by an isotropic distribution of velocities 
and a Nuker (Lauer et al 1995) density profile p(r) 
with specified parameter values. 

• Insert a black hole binary with specified masses 
Mi = M 2 = M and radius r/,. For 'realistic' values 
of M and r^, M(r/,) is typically small compared with 
the black hole mass, so that one can assume, at least 
approximately, that the black holes are executing a 
Keplerian orbit with frequency uj = ^/ Mj 4r^ . 

• Next evolve the initial conditions in the fixed time- 
dependent potential comprised of the Nuker potential 
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plus the potential of the orbiting black hole binary, 
and track the radial density distribution p(r,t). 
• Finally, assuming that mass traces light, compute 
line-of-sight integrals along the density distribution to 
obtain integrated surface densities and, hence, surface 
brightness distributions as functions of time. 

Although this approach does not pretend to be 
completely realistic, it would not seem totally unrea- 
sonable to insert the binary 'by hand' without allow- 
ing for the dynamics whereby it has evolved into a 
tightly bound orbit near the galactic center. When 
the binary orbit is very large, it will have a compar- 
atively minimal effect. Energy and mass transport 
only becomes important at comparatively small radii, 
where M > M(r/»), and again become unimportant 
when the radius becomes too small. Most of the ac- 
tion happens for a relativley limited range of radii. 

Note, moreover, that the assumption M J> M{rh) 
tends to mitigate the fact that the computations are 
not fully self-consistent: Although the bulk forces 
associated with the galaxy cannot be neglected at 
all radii where the binary has an appreciable effect, 
they can presumably be neglected, at least approxi- 
mately, at the comparatively small radii near the bi- 
nary where the effect of the black holes is strongest. 

4.2. The initial form of the density and po- 
tential 

Initial attempts at modeling using a spherical Dehnen 
potential yielded results in qualitative agreement with 
observations. However, comparatively large system- 
atic deviations were observed, which appeared to re- 
flect the fact that the transition between the inner 
and outer power-law profiles predicted by a Dehnen 
potential is too gradual to represent real galaxies. For 
this reason, models were constructed instead using an 
initial density distribution satisfying the more general 
Nuker Law (Lauer et al 1995) 



Po(r) = p c r^(l + r a ) 



(10) 



Dehnen models are recovered for a = 1 and (3 = 4. 
The central density p c was chosen so that the total 
galactic mass M g = 1.0. The associated potential 
V(r) satisfies (in units with G = 1) 

- / p(r)r 2 df + / p(f)fdf 
r Jo Jr 



V(r) = -4tt 



(11) 



Unfortunately, this potential can be expressed in 
terms of elementary functions only for certain choices 



of a and j3, which means, generically, that orbits 
must be computed using an expensive interpolation 
scheme. This motivated an effort to seek fits assum- 
ing values of a and (3 for which V can be expressed 
analytically. For the small number of profiles consid- 
ered hitherto, reasonable fits were achieved for a = 2 
and (3 = 4, which, for 7 — 0, yields a potential 



V(r) = - 



2 tan 



and 



M(r 



tan 



(12) 



(13) 



Most models were constructed assuming M(r^) <C 
Mi + M2, so that the approximation of a Kcplerian 
frequency is typically very good. However, in an effort 
to allow for the influence of the galactic potential, 
the models allowed for a slightly modified frequency 
lu = y/M/(2r h ) 3 , where M = M 1 +M 2 + 4M(r h ). 

4.3. Generating a surface brightness distribu- 
tion 

Configuration space was divided into N = 100 
equally spaced concentric shells i. Each shell cor- 
responded to a range of energies, I^-i < E < Ei, 
i = 1,...,N, sampled along the principal axes in the 
plane of the binary, but perpendicular to the line con- 
necting them. This was done to ensure that energy 
was a monotonic function of radius, so that shuffling 
of energies could be related directly to a redistribu- 
tion of orbits in configuration space. Each shell was 
sampled to select M — 300 initial conditions, which 
were integrated for a time t = 512. Orbital data were 
recorded periodically and the new energies used to 
reassign orbits to (in general) new shells. If Mi(t) 
denotes the number of orbits in shell i at time t, then 



Ai(t) = 



M ' 



(14) 



the relative fluctuation in number, can be interpreted 
as a discretised version of a radial density fluctuation 
5(t) satisfying 



p{r,t) = [l + 8{r,t)]p (r), 



(15) 



with po the initial density. S(t) was interpolated from 
Ai(t) using a smooth-curve fitting routine. 

The resulting density p(r, t) was then integrated 
along the line of sight to generate a surface brightness 



p(r,t) 



2 r 



p(f, t)f 

\Jf 2 — r 2 



dr. 



(16) 
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Hero T denotes the mass-to-light ratio, which was 
assumed constant for the modeling described here. 

4.4. Results 

Figure 13 exhibits data for a typical model, cor- 
responding to a Nuker Law with a = 2, /3 = 4, 
and 7 = 0. The binary parameters are M = 0.005, 
r/j = 0.15, and to = 0.6086. The half-mass radius is 
r = 2.264; 75% of the mass is contained with r = 5. 

It is obvious that the binary induces a distinctive 
signature, characterised by an inversion in both the 
mass density and the surface brightness profile. In 
some cases, especially when 7^0, the contents of the 
innermost shells can remain essentially intact. Aside, 
however, from those innermost shells, one can identify 
a well-defined sphere of influence where the binary has 
observable effects. For r < r\, there is a systematic 
undcrpopulation of stars and, hence, a dip in lumi- 
nosity; for n < r < r 2 , there is a systematic over- 
population or bulge resulting from stars transported 
outwards from radii r < r 2 . For r > r 2 the density 
and surface brightness distribution remain essentially 
unchanged. Significantly, the signature, once estab- 
lished, remains largely unchanged in bulk properties. 
In particular, the dip is comparably prominent visu- 
ally at t = 128 and t = 512. 

Figure 14 exhibits a more systematic attempt to 
model the surface brightness of NGC 3706, again 
starting from a Nuker Law with a = 2, (3 = 4, and 
7 = 0. Physical distance was translated into angu- 
lar separation assuming a scaling such that r = 1 
corresponds to 0.15 arcsec or, given the distance es- 
timate given by Lauer et al, r 24 pc. Once again 
M = 0.005, but now rt = 0.085, which corresponds 
to a physical ~ 2.0 pc and an angular separation 
~ 0.014 arcsec. The orbital frequency u> — 1.567, 
which implies an orbital period r w 4.00. 

In this case, the inner dip extends out to <~ 0.10 
arcsec; the bulge extends to <~ 0.30 arcsec. On scales 
;> 0.30 arcsec, i.e., r ^ 75 pc, the binary has only 
a comparatively minimal effect, so that the surface 
brightness remains essentially unchanged. 

The perturbed Nuker model is quite successful in 
modeling the dip and the outer region, where errors in 
surface brightness correspond typically to 5/i ~ 0.005 
magnitudes or less. In particular, the dip is much 
better fit by the perturbed model than by an unper- 
turbed Nuker model. Both qualitatively - in that an 
unperturbed Nuker model requires a monotonically 



decreasing surface brightness - and quantitatively - 
in terms of the actual error S/j, -, the perturbed model 
does a much better job. However, both the perturbed 
and unperturbed models are somewhat less success- 
ful in accounting for the detailed shape of the bulge 
(although the model with a binary does a bit better). 

There are at least two possible explanations for 
this lack of success. Most obvious is the fact that, 
demanding a = 2 and (i = 4, so that the poten- 
tial could be written in terms of elementary func- 
tions, limits one's flexibility in modeling the transition 
region between the inner and outer (unperturbed) 
power law profiles. Allowing for fractional parameter 
values (which requires that the potential be computed 
numerically) will likely yield better fits. However, it 
is also possible that this lack of success reflects in part 
the oversimplistic character of this kinematic model. 
In a real galaxy, the binary orbit decays as the bi- 
nary tranfers energy to the stars; and the fact that 
Th is not really constant might be expected to have 
some observable effects. Attempts to remedy these 
deficiencies of the model are currently underway. 

It should, however, be stressed that the general 
effects of the binary are relatively insensitive to r/j, 
provided only that M > M(rh)- This is, e.g., evident 
from Fig. 15, which exhibits surface brightness dis- 
tributions at t = 256 for both the model considered 
in Fig. 14 and another model identical except that 
r/j = 0.025, a radius only 0.29 times as large. There 
are some differences in detail, but neither fit is clearly 
superior visually. It is also evident from Fig. 14 that 
the basic observable structure develops very quickly. 
Although the details of the surface brightness profile 
can vary considerably for times as long as t ~ 128 or 
more, the existence of the dip region is obvious al- 
ready by t ~ 32, about 8 binary orbital periods for 
the r h = 0.085 model. 

Although the model described here is kinematic, 
one can try to describe how it might be manifested 
in a self-consistent description: When the binary or- 
bit is too large, its decay will be dominated by more 
conventional processes, discussed, e.g., in Tremaine 
& Weinberg (1984) and Nelson & Tremaine (1999). 
However, once the radius is sufficiently small that 
M(r/j) ~ 2M, the resonant phase mixing described 
here - which can be viewed as a variant of Tremaine's 
resonant relaxation - will be triggered. As additional 
energy is transferred to the stars, the binary will con- 
tinue to decay and, when the size of the orbit becomes 
too small, the effect will again 'turn off.' 
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For the models in Fig. 14 and 15, the process 
should 'turn on' at Th,\ ~ 0.2 and 'turn off' at r/,,2 ~ 
0.005. However, during this interval, the binary will 
have lost an energy ~ M 2 /rh.2 ~ 0.02, several per- 
cent of the energy of the galaxy at the time that the 
process begins. The obvious questions, therefore, are: 
How long does it take for a binary with radii satis- 
fying rh.2 <z tk <, r h,i to pump this much energy into 
the stars? And is this long enough to establish the 
signature observed in Figs. 14 and 15? The time re- 
quired depends to a certain extent on the precise value 
of Th- However, an analysis of the models in Figs. 14 
and 15, as well as models with somewhat larger and 
smaller values of rt indicates that the total energy 
required to establish the observed signature is rela- 
tively small. For example, the model with rh = 0.085 
entailed an increase in galactic energy of order 1% at 
t = 32 and 3% at t = 512. The model with r h = 0.025 
yielded 1.5% at t = 32 and 6% at t = 512. 

Alternatively, a binary decay rate can be estimated 
as follows: Given that the pumping of energies into 
the stars is diffusive, the decay time r should sat- 
isfy t/T ~ (M 2 /r h (SE)) 2 , where T is the time over 
which (SE) is computed. Supposing, however, that 
r h ~ 0.01, M(r h ) ~ 0.01, (SE) ~ 0.01, and T ~ 10, 
one infers that t ~ T ~ 10. For the orbit to shrink 
from rh = 0.2 to r/, = 0.005, a factor of 40, would 
require a few r, say t ~ 50, an interval long enough 
to establish a distinctive luminosity dip. 

One final point should be noted. Attributing such a 
luminosity dip to a supermassive binary does not nec- 
essarily imply that the binary should still be present. 
If, neglecting the binary, the galaxy can be idealised 
as a collisionless equilibrium, one might expect that a 
dip in surface brightness, once generated, could per- 
sist even after the binary has coalesced, at least for 
times short compared with the time scale on which 
stars at larger radii could be scattered inwards via 
collisional relaxation. To the extent that the bulk po- 
tential is time-independent, in the absence of 'colli- 
sions' energy is conserved, so that an underpopulated 
region in energy space cannot be repopulated. 

5. Discussion 

The computations described here yield several sig- 
nificant conclusions about phase mixing in a time- 
dependent potential. Most obvious is the fact that a 
supermassive black hole binary can serve as an im- 
portant source of transient chaos which facilitates ef- 



ficient resonant phase mixing, shuffling the energies of 
stars (or any other objects) as well as phase space co- 
ordinates on a constant energy hypersurface. In par- 
ticular, the effects observed here from a comparatively 
'small scale' perturbation are quite similar to the ef- 
fects observed when galaxies are subjected to larger 
scale systematic oscillations (Kandrup, Vass & Sideris 
2003). It is especially striking that, even though the 
perturbation is relatively low amplitude and concen- 
trated very near the center of the galaxy, it can have 
significant effects at comparatively large radii. All 
this reinforces the expectation that resonant phase 
mixing could be a generic physical effect in galaxies 
subjected to an oscillatory time dependence. 

Contrary, perhaps, to naive expectation, it appears 
that the shuffling of energies is diffusive, rather than 
exponential, so that energy phase mixing is less dra- 
matic than phase mixing of coordinates and velocities. 
Even though the time-dependent perturbation can in- 
crease both the relative abundance of chaotic orbits 
and the degree of exponential sensitivity exhibited by 
chaotic orbits, it need not make orbits exponentially 
unstable in the phase space direction orthogonal to 
the constant energy hypersurfaces. 

However, such energy shuffling could still play an 
important role in violent relaxation. Indeed, the fact 
that this energy shuffling is not exponential is con- 
sistent with self-consistent simulations of violent re- 
laxation (e.g., Quinn & Zurek 1988) which indicate 
that, even though 'particles' are almost completely 
'randomised' in terms of most phase space coordi- 
nates, they exhibit a partial remembrance of initial 
conditions. In particular, 'particles' that start with 
low (high) binding energies tend systematically to end 
with low (high) binding energies. If, e.g., stars in sim- 
ulations involving hard, head-on collisions of galaxies 
are ordered in terms of their initial and final binding 
energies, the rank correlation 1Z between the initial 
and final ordered lists typically satisfies (Kandrup, 
Mahon & Smith 1993) K £ 0.6. 

That a supermassive binary will cause a system- 
atic readjustment in the density distribution of the 
host galaxy seems largely independent of the form 
of the galactic potential or the orbital parameters of 
the binary, although the precise form of the readjust- 
ment does depend on these details. In particular, one 
sees qualitatively similar effects for Dehnen potentials 
with different cusp indices 7 and for Nuker Laws with 
different transitional radii properties. Similarly, the 
eccentricity and the orientation of the supermassive 
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binary are not crucial, and allowing for unequal, but 
still comparable, masses does not result in qualitative 
changes. Irrespective of all these details, when the to- 
tal binary mass Mi + Mi <C M{rh), with rt the 'size' 
of the binary orbit, stars cannot resonate with the bi- 
nary and comparatively little mass transport occurs. 
However, when Mi + Mi ~ M(rh), one starts seeing 
substantial effects which can extend to radii » tv 

One might therefore expect that, when its orbit is 
large, the binary will have only a minimal effect on 
the bulk properties of the galaxy; but that when, as 
a result of dynamical friction (e.g., Merritt 2001), the 
orbit has decayed to a sufficiently small size, it will 
begin to have an appreciable - and observable - effect. 
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Fig. 1. — (a) The mean shift in energy (SE) for all 
the orbits in a 1600 orbit ensemble with E — 0.87 and 
(fin) ~ 0.86, evolved in a spherical oscillator potential 
with M = 0.05, r h = 0.3, and a 2 = b 2 = c 2 = 1.0, 
plotted as a function of frequency u>. (b). The dis- 
persion (TgE for all the orbits, (c) - (d) The same as 
the preceding for orbits integrated in a potential with 
a 2 = 1.33, b 2 = 1.0, and c 2 = 0.80 and an ensemble 
with E = 0.87 and (r rn ) w 0.89. 
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Fig. 2. — (a) The mean energy shift (SE) for the same 
ensemble used to generate FIG. 1 (c) and (d), inte- 
grated with r h = 0.3, a 2 = 1.33, b 2 = 1.0, c 2 = 0.80, 
and (from top to bottom) M = 0.05, M = 0.0281, 
M = 0.0158, and M = 0.005. (b) (SE) for the same 
ensembles - solid line for M — 0.005, dashed line for 
M = 0.0158, dot-dashed line for M = 0.0281, and 
triple-dot-dashed for M = 0.05 - now expressed in 
units of the maximum shift {5E max ) . (c) The mean 
energy shift (SE) for integrations with a 2 = 1.33, 
b 2 = 1.0, c 2 = 0.80, M = 0.05, and (curves peak- 
ing from left to right) r/, = 0.4. r/, = 0.3, = 0.2, 
and rh = 0.1. (d) (<5i?) expressed in units of the max- 
imum shift (SE max ) - solid line for r/, = 0.1, dashed 
for r/j = 0.2, dot-dashed for = 0.3, and triple-dot- 
dashed for rh = 0.4. 
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Fig. 3. — (a) The fraction / of strongly chaotic or- 
bits in a 1600 orbit ensemble with initial energy E = 
—0.70 and mean initial radius (rj n ) w 0.33, evolved in 
a spherically symmetric Dehnen potential with 7 = 
1.0, M = 0.01, r = 0.05, and a 2 = b 2 = c 2 = 1.00. 
(b) The mean value (x) of the largest finite time Lya- 
punov exponent for the strongly chaotic orbits, (c) 
The mean shift in energy (SE) for all the orbits, (d) 
The dispersion <tse) for all the orbits, (e) - (h) The 
same as the preceding for orbits integrated in a poten- 
tial with a 2 = 1.25, b 2 = 1.00, and c 2 = 0.75, again 
for an ensemble with E = —0.70 and (r in ) m 0.33. 
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Fig. 4. — (a) a 2 E , where asE(t) is the time-dependent 
spread in energy shifts associated with an ensemble 
of orbits evolved in an oscillator potential with M = 
0.05, r h = 0.3, a 2 = 1.33, b 2 = 1.0, c 2 = 0.80 and 
u) = 0.5. (b) (SE(t)) 2 for the same ensemble, (c) and 
(d) The same for w = 1.0. (e) and (f) The same for 
uo = 2.0. (g) and (h) The same for uj — 4.0. (i) and 
(j) The same for u — 8.0. 
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Fig. 5. — (a) Scatter plots relating ose and \-> where 
a be represents the dispersion associated with the 
time-dependent 8E{t) for an individual orbit over the 
interval < t < 512. The orbits are the same that 
were used to generate FIG. 3, integrated with to = 0.5. 
(b) Scatter plots relating (SE) and x, where (SE) rep- 
resents the mean value of SE(t), computed for the 
same orbits as in (a), (c) and (d) The same for 
lu = 1.0. (c) and (f) The same for lu = 2.0. (g) 
and (h) The same for ui = 4.0. (i) and (j) The same 
for u = 8.0. 



0.4 

0.2 

>- 0.0 

-0.2 
-0.4 



(o) 



0.4 

0.2 

>- o.o 

-0.2 
-0.4 



-0.4 -0.2 0.0 0.2 0.4 

X 



-0.4 -0.2 0.0 0.2 0.4 

X 



0.-1 




0.4 


0.2 




0.2 


0.0 




>- 0.0 


0.2 




-0.2 


0.4 


(b) 


-0.4 



(g) 



-0.4 -0.2 0.0 0.2 0.4 
X 




0.2 

o. 'c; • 

-0.4 -0.2 0.0 0.2 0.4 



0.4 

0.2 
>- 0.0 
-0.2 



-0.4 L 

-0.4 -0.2 0.0 0.2 0.4 



0.4 

0.2 
0.0 
-0.2 
-0.4 



(j) 



-0.4 -0.2 0.0 0.2 0.4 

X 



Fig. 6. — (a) The x and y coordinates at t = for an 
initially localised ensemble of orbits with E = —0.70 
and (ri n ) « 0.33, evolved in a spherical Dehnen po- 
tential with 7 = 1.0, r h = 0.05, M = 0.005, and 
lu = VlO. (b) t = 8. (c) t = 16. (d) t = 32. (e) 
t = 64. (f) - (j). The same for stationary black holes, 
i.e., to = 0.0. 
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Fig. 7. — (a) The mean shift in energy, (SE), com- 
puted for an ensemble of orbits with E = —0.70 and 
(r) w 0.33, evolved in a 7 = 1.0 Dehnen model with 
a 2 = b 2 = c 2 = 1 in the presence of a supermas- 
sive binary comprised of two black holes executing 
strictly circular orbits with ru = 0.05 and different 
values of Mi = M 2 = M. (b) (SE) for the same 
ensemble evolved in the same Dehnen model, again 
allowing for a binary executing circular orbits, but 
now with Mi = M 2 = 0.01 and variable rv (c) and 
(d) The same for a model with a 2 = b 2 = 0.90 and 
c 2 = 1.21. (e) and (f) The same for a model with 
a 2 = b 2 = 1.21 and c 2 = 0.64 (g) and (h) The same 
for a model with a 2 = 1.10, b 2 = 1.0 and c 2 = 0.90. 
(i) and (j) The same for a model with a 2 = 1.25, 
b 2 = 1.0 and c 2 = 0.75. In each case, the frequency 
lj = yj (Mi + M 2 )/a 3 , with A the semi-major axis. 
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Fig. 8. — (a) The mean shift in energy, (SE), com- 
puted for an ensemble of orbits with E = —0.70 and 
(r) w 0.33, evolved in a 7 = 1.0 Dehnen model with 
a 2 = b 2 = c 2 = 1 in the presence of a supermas- 
sive binary comprised of two black holes with mass 
Mi = M 2 = 0.01 executing orbits with semi-major 
axis A = 0.10 and variable eccentricity e. (b) (SE) 
for the same ensemble evolved in the same Dehnen 
model, again assuming M\ + M 2 = 0.02 and a = 0.10, 
but now allowing for different ratios M 2 j(M\ + M 2 ). 
(c) and (d) The same for a model with a 2 = b 2 = 0.90 
and c 2 = 1.21. (e) and (f) The same for a model with 
a 2 = b 2 = 1.21 and c 2 = 0.64 (g) and (h) The same for 
a model with a 2 = 1.10, b 2 = 1.0 and c 2 = 0.90. (i) 
and (j) The same for a model with a 2 = 1.25, b 2 = 1.0 
and c 2 = 0.75. 
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Fig. 9. — (a) The initial angle-averaged radial density 
distribution associated with a 4800 orbit sampling of 
the E = —0.70 constant energy hypersurface, subse- 
quently integrated in a pseudo-Dehnen potential with 
7 = 1.0, M = 0.01, r h = 0.05, a 2 = 1.25, b 2 = 1.00, 
c 2 = 0.75 and uj = V^O. (b) The density at t = 16. 
(The dotted line reproduces the initial distribution.) 
(c) t = 32. (d) t = 64. (e) t = 128. (f) - (j) The same 
for stationary black holes, i.e., u = 0.0. 
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Fig. 10. — (a) Mean energy shift (5E), computed for 
ensembles with different radii, for orbits in a spherical 
Dehnen model with 7 = 0.0 and a = b = c = 1.0 and 
black hole parameters M = 0.005, = 0.25, and 
uo = 0.2828. Note that the radius r is plotted on a 
logarithmic scale, (b) The dispersion <jse for the same 
ensembles, (c) The mean value (x) for each ensemble, 
(d) The dispersion <r x . (e) - (h) The same for a model 
with 7 = 1.0. 
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Fig. 11. — (a) The direction-dependent spatial dis- 
tributions (solid curve) and ro(|z|) (dashes) at 
t = 512 generated for a 4800 orbit sampling of the 
E = -0.70 hypcrsurfacc with 7 = 1.0, M = 0.01, 
r/j = 0.005, a = b = c = 1, and uo = \/20, along with 
the distribution n(|x|) (dot-dashed) at time t = 0. (b) 
The corresponding direction-dependent velocity dis- 
tributions, (c) x and y coordinates for the ensemble 
at t = 512. (d) x and z coordinates for the ensemble 
at t = 512. 
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Fig. 12. — Direction-dependent spatial distributions 
n(|a;|) (solid curve), n(|y|) (dashes), and n(|z|) (dot- 
dashes) generated for a 4800 orbit sampling of the 
E = -0.62 hypersurface with 7 = 1.0, M = 0.01, 
r h = 0.005, a 2 = 1.25, b 2 = 1.0, c 2 = 0.75, and 
uo = V20. along with the distribution n{\x\) (dot- 
dashed) at time t — 0. (a) The distributions at time 
t = 0. (b) The distributions at £ = 512, allowing for 
a binary orbiting in the x — y plane, (c) The same for 
a binary in the y — z plane, (d) The z — x plane. 
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Fig. 13. — Computed quantities for a Nuker model 
with a = 2, (3 = 4, 7 = 0, Mbh = 0.005, r h = 
0.15, and u> = 0.6086. The first column exhibits A, 
the relative fluctuation in density for different shells; 
the second exhibits the interpolated smooth density 
p; the third exhibits the surface brightness, assuming 
that mass traces light. From top to bottom, rows 
represent integration times t = 128, 256, 384, and 512. 
In each case, the dotted lines represent the original 
unperturbed values. 




Fig. 14.— Modeling NGC 3706 with a Nuker model 
with a = 2, (3 = 4, 7 = 0, r h = 0.085, and cu = 1.567. 
The left column exhibits the observed surface bright- 
ness profile (solid circles), the surface density pre- 
dicted by an unperturbed Nuker Law (dotted lines), 
and the time-dependent surface density generated by 
the binary (solid lines) at times (from top to bottom) 
t = 32, t = 64, t = 128, and t = 256. The right 
column exhibits the relative error of the the fit for a 
Nuker model without (dashes) and with the binary 
(solid lines). 
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Fig. 15.— (a) The best fit model with a = 2, (3 = 4, 
7 = 0, r h = 0.025, and uj = 8.968 at time t = 256. (b) 
The same, except assuming — 0.085 and uu — 1.567. 
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